#install.packages('ggtree')
#install.packages('ggplot2')
#install.packages('dplyr')
#install.packages('treeio')
#install.packages('phytools')
#install.packages('ape')
#install.packages('ggpubr')

library(ggtree)
## Registered S3 method overwritten by 'treeio':
##   method     from
##   root.phylo ape
## ggtree v2.2.1  For help: https://yulab-smu.github.io/treedata-book/
## 
## If you use ggtree in published research, please cite the most appropriate paper(s):
## 
## - Guangchuang Yu. Using ggtree to visualize data on tree-like structures. Current Protocols in Bioinformatics, 2020, 69:e96. doi:10.1002/cpbi.96
## - Guangchuang Yu, Tommy Tsan-Yuk Lam, Huachen Zhu, Yi Guan. Two methods for mapping and visualizing associated data on phylogeny using ggtree. Molecular Biology and Evolution 2018, 35(12):3041-3043. doi:10.1093/molbev/msy194
## - Guangchuang Yu, David Smith, Huachen Zhu, Yi Guan, Tommy Tsan-Yuk Lam. ggtree: an R package for visualization and annotation of phylogenetic trees with their covariates and other associated data. Methods in Ecology and Evolution 2017, 8(1):28-36. doi:10.1111/2041-210X.12628
library(ggplot2)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(treeio)
## treeio v1.12.0  For help: https://yulab-smu.github.io/treedata-book/
## 
## If you use treeio in published research, please cite:
## 
## LG Wang, TTY Lam, S Xu, Z Dai, L Zhou, T Feng, P Guo, CW Dunn, BR Jones, T Bradley, H Zhu, Y Guan, Y Jiang, G Yu. treeio: an R package for phylogenetic tree input and output with richly annotated and associated data. Molecular Biology and Evolution 2020, 37(2):599-603. doi: 10.1093/molbev/msz240
library(phytools)
## Loading required package: ape
## 
## Attaching package: 'ape'
## The following object is masked from 'package:treeio':
## 
##     drop.tip
## The following object is masked from 'package:ggtree':
## 
##     rotate
## Loading required package: maps
## 
## Attaching package: 'phytools'
## The following object is masked from 'package:treeio':
## 
##     read.newick
library(ape)
library(ggpubr)
## 
## Attaching package: 'ggpubr'
## The following object is masked from 'package:ape':
## 
##     rotate
## The following object is masked from 'package:ggtree':
## 
##     rotate
setwd('/Users/owlex/Dropbox/Documents/Northwestern/rcs_consult/r_phylogenetics_workshop/r_phylogenetics_worshop') # change this to your R Markdown's file path

We will learn about reading and plotting trees using a set of 20 genes that belong to four different classes of enoyl reductases (five per class) that are involved in fatty acid synthesis. These are fabK, fabG, fabV, and fabI.

For your own interpretation, know that FabV, FabK, and FabI have the same enzymatic function, but fabG does not. Despite the overlap in function, fabG and FabI are more closely related to each other by sequence. FabV is the next most closely related, and fabK is the least related.

Throughout the workshop, I highlight when I’m using non-ggtree packages by calling on the package name first and then the function. For example using the ape package’s function read, I call on it by typing ape::read(). Hopefully this makes the roles of the different packages being used more obvious and how they fit together more clear.

Section 1: Reading in different phylogenetic tree file types

There are a variety of different phylogenetic tree file types. One of the most popular ones is the .newick (also abbreviated .nwk)

We can read in the file using ape’s read.newick() function.

tree_nwk <- ape::read.tree('raw_enoyl_seqs.nwk')

tree_nwk is of class ‘phylo’. When can always confirm this by using the class() function.

class(tree_nwk)
## [1] "phylo"

We can access the information about the tree object

# general info
tree_nwk
## 
## Phylogenetic tree with 20 tips and 18 internal nodes.
## 
## Tip labels:
##  tr|B0BP34|B0BP34_ACTPJ, tr|A0A1R0ED53|A0A1R0ED53_HAEPA, tr|A0A0M3G239|A0A0M3G239_HAEHA, sp|P0AEK4|FABI_ECOLI, sp|O24990|FABI_HELPY, sp|I6Y778|FABG4_MYCTU, ...
## Node labels:
##  , 0.197, 0.679, 0.990, 0.995, 0.183, ...
## 
## Unrooted; includes branch lengths.
#  tip label names
tree_nwk$tip.label
##  [1] "tr|B0BP34|B0BP34_ACTPJ"         "tr|A0A1R0ED53|A0A1R0ED53_HAEPA"
##  [3] "tr|A0A0M3G239|A0A0M3G239_HAEHA" "sp|P0AEK4|FABI_ECOLI"          
##  [5] "sp|O24990|FABI_HELPY"           "sp|I6Y778|FABG4_MYCTU"         
##  [7] "sp|P0AEK2|FABG_ECOLI"           "sp|P9WGT3|MABA_MYCTU"          
##  [9] "sp|P71534|MABA_MYCS2"           "sp|P24182|ACCC_ECOLI"          
## [11] "sp|Q73Q47|FABV_TREDE"           "sp|Q9HZP8|FABV_PSEAE"          
## [13] "sp|Q62L02|FABV_BURMA"           "sp|Q9KRA3|FABV1_VIBCH"         
## [15] "sp|Q2P9J6|FABV_XANOM"           "sp|P54616|FABI_BACSU"          
## [17] "sp|Q6GI75|FABI_STAAR"           "sp|Q81GI3|FABI_BACCR"          
## [19] "tr|E6KYQ2|E6KYQ2_9PAST"         "tr|A0A4Y9GUQ1|A0A4Y9GUQ1_9PAST"

The ggtree package makes it very easy to quickly read in a .newick file and view a nicely formatted phylogenetic tree.

tree_nwk <- ape::read.tree('raw_enoyl_seqs.nwk')
ggtree(tree_nwk)+ # read in tree
  geom_tiplab()+ # show tip labels
  geom_treescale(5) # recenter the tree
## Warning: `tbl_df()` is deprecated as of dplyr 1.0.0.
## Please use `tibble::as_tibble()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.

Not all phylogenetic tree files can be read with read with ape::read.tree(). The treeio package can open virtually all other existing tree file types.. Click on the provided link see their commands.

Section 2: Basic tree plotting with ggtree

Another common tree filetype are those produced by the program RAxML. These are a modified newick file that contains bootstrap numbers at each node.

We will use the ape::read.tree() to open a raxml file in R.

Exercise 1

  1. Use ape::read.tree() to read in the file raw_enoyl_seqs.nwk and assign it to the variable my_nwk

  2. Display the my_nwk tree using ggtree()

  3. Show the geom_tiplab() for each leaf in the ggtree() visualization of my_nwk

  4. Modify the code from part 3 so that the x-axis is rescaled at x=7 using geom_treescale()

Some file types hold additional information. For example, the RAxML tree files can contain bootstrap values at each node bifurcation.

tree_raxml <- ape::read.tree('RAxML_bipartitions.raw_enoyl_seqs')
class(tree_raxml)
## [1] "phylo"
# for a RAxML file, the bootstrap values are stored under the node.label 
tree_raxml$node.label
##  [1] ""    "100" "54"  "100" "56"  "94"  "47"  "99"  "77"  "97"  "100" "48" 
## [13] "29"  "46"  "100" "85"  "100" "55"

The bootstrap values can be displayed “out of the box” on the tree using the option geom_nodelabaes(label=label)

ggtree(tree_raxml)+
  geom_tiplab()+
  geom_nodelab(aes(label=label))+ #show bootstraps
  geom_treescale(x=5)

Oftentimes we want to re-order the tree according to a phylogenetic rooting.

The phytools package provides a midpoint.root()ing function for re-ordering the tip labels so that the root is at the two longest branchest.

tree_raxml_midpointroot <- phytools::midpoint.root(tree_raxml)

We can compare the order of the original tree vs the midpoint by accessing data stored in the tree object

compare_tips <- data.frame('original'=tree_raxml$tip.label, 
                           'midpoint'=tree_raxml_midpointroot$tip.label)

compare_tips
##                          original                       midpoint
## 1            sp|Q9HZP8|FABV_PSEAE           sp|O24990|FABI_HELPY
## 2            sp|P24182|ACCC_ECOLI           sp|P54616|FABI_BACSU
## 3            sp|O24990|FABI_HELPY           sp|Q81GI3|FABI_BACCR
## 4            sp|P54616|FABI_BACSU           sp|Q6GI75|FABI_STAAR
## 5            sp|Q81GI3|FABI_BACCR           sp|P0AEK4|FABI_ECOLI
## 6            sp|Q6GI75|FABI_STAAR tr|A0A0M3G239|A0A0M3G239_HAEHA
## 7            sp|P0AEK4|FABI_ECOLI tr|A0A1R0ED53|A0A1R0ED53_HAEPA
## 8  tr|A0A0M3G239|A0A0M3G239_HAEHA tr|A0A4Y9GUQ1|A0A4Y9GUQ1_9PAST
## 9  tr|A0A1R0ED53|A0A1R0ED53_HAEPA         tr|E6KYQ2|E6KYQ2_9PAST
## 10 tr|A0A4Y9GUQ1|A0A4Y9GUQ1_9PAST         tr|B0BP34|B0BP34_ACTPJ
## 11         tr|E6KYQ2|E6KYQ2_9PAST          sp|I6Y778|FABG4_MYCTU
## 12         tr|B0BP34|B0BP34_ACTPJ           sp|P0AEK2|FABG_ECOLI
## 13          sp|I6Y778|FABG4_MYCTU           sp|P71534|MABA_MYCS2
## 14           sp|P0AEK2|FABG_ECOLI           sp|P9WGT3|MABA_MYCTU
## 15           sp|P71534|MABA_MYCS2           sp|Q9HZP8|FABV_PSEAE
## 16           sp|P9WGT3|MABA_MYCTU           sp|P24182|ACCC_ECOLI
## 17           sp|Q73Q47|FABV_TREDE           sp|Q73Q47|FABV_TREDE
## 18           sp|Q2P9J6|FABV_XANOM           sp|Q2P9J6|FABV_XANOM
## 19          sp|Q9KRA3|FABV1_VIBCH          sp|Q9KRA3|FABV1_VIBCH
## 20           sp|Q62L02|FABV_BURMA           sp|Q62L02|FABV_BURMA

We can also compare them by plotting both trees side by side

# tree 1

tree_original_visual <- ggtree(tree_raxml)+
  geom_tiplab()+
  geom_nodelab(aes(label=label))+
  geom_treescale(x=10)
  

tree_midpoint_visual <- ggtree(tree_raxml_midpointroot)+
  geom_tiplab()+
  geom_nodelab(aes(label=label))+
  geom_treescale(x=10)
  

ggpubr::ggarrange(tree_original_visual, tree_midpoint_visual, nrow=1)

The tree on the right is far easier to understand. The different gene types almost cluster together in monophyletic branches.

Another option to re-organize the tree is to root the tree to an ancestral gene or to an outlier.

We will try this by using ape::root() on the gene sp|P24182|ACCC_ECOLI and tr|E6KYQ2|E6KYQ2_9PAST

# first outgroup tree
tree_raxml_outgrouproot <- ape::root(tree_raxml ,outgroup='sp|P24182|ACCC_ECOLI')

tree_outgroup_visual <- ggtree(tree_raxml_outgrouproot)+
  geom_tiplab()+
  geom_nodelab(aes(label=label))+
  geom_treescale(x=10)


# second outgroup tree
tree_raxml_outgrouproot2 <- ape::root(tree_raxml ,outgroup='tr|E6KYQ2|E6KYQ2_9PAST')

tree_outgroup2_visual <- ggtree(tree_raxml_outgrouproot2)+
  geom_tiplab()+
  geom_nodelab(aes(label=label))+
  geom_treescale(x=10)

# plot together
ggpubr::ggarrange(tree_outgroup_visual,tree_outgroup2_visual,nrow=1)

The second tree provides a similar ordering that accentuates the divergence of fabV.

Exercise 2

  1. Use ape::read.tree() to read in the tree in 'RAxML_bipartitions.raw_enoyl_seqs' and store it to the variable my_tree_raxml.

  2. Use ape::root() to root the my_tree_raxml to the outgroup = 'sp|P0AEK4|FABI_ECOLI' and store it in the variable my_rooted_tree

  3. Visualize the my_tree_raxml using ggtree() with with annotation layers geom_treescale(x=10), geom_tiplab(), and geom_nodelab(aes(label=label))

The ggtree package offers a lot of customization options for your tree plots similar to ggplot2.

We will now go over a variety of them below.

First we will make a tree for a base comparison

tree_raxml <- ape::read.tree('RAxML_bipartitions.raw_enoyl_seqs')

tree_raxml_outgrouproot_compare <- ape::root(tree_raxml ,outgroup='tr|E6KYQ2|E6KYQ2_9PAST')

tree_raxml_outgrouproot_visual <- ggtree(tree_raxml_outgrouproot_compare)+
  geom_tiplab()+
  geom_nodelab(aes(label=label))+
  geom_treescale(x=5)

Making tip labels smaller, making bootstrap values smaller and moving them a bit so they over lap less

tree_compare <- ggtree(tree_raxml_outgrouproot_compare)+
  geom_tiplab(size=2)+
  geom_nodelab(aes(label=label),nudge_x=-.18,nudge_y=.4,size=2)+
  geom_treescale(x=5)

ggpubr::ggarrange(tree_raxml_outgrouproot_visual,tree_compare,nrow=1)

Adding tip points and offsetting the tip labels to accomodate them

tree_compare <- ggtree(tree_raxml_outgrouproot_compare)+
  geom_tiplab(size=2,offset=0.1)+
  geom_tippoint()+
  geom_nodelab(aes(label=label),nudge_x=-.18,nudge_y=.4,size=2)+
  geom_treescale(x=7)

ggpubr::ggarrange(tree_raxml_outgrouproot_visual,tree_compare,nrow=1)

Convert the tree into a cladogram

tree_compare <- ggtree(tree_raxml_outgrouproot_compare, branch.length = 'none')+
  geom_tiplab(size=2,offset=0.3)+
  geom_tippoint()+
  geom_nodelab(aes(label=label),nudge_x=-.6,nudge_y=.4,size=2)+
  geom_treescale(x=25)

ggpubr::ggarrange(tree_raxml_outgrouproot_visual,tree_compare,nrow=1)

Have all tip labels be right-justified

tree_compare <- ggtree(tree_raxml_outgrouproot_compare)+
  geom_tiplab(size=2,align = TRUE,offset=.3)+
  geom_tippoint()+
  geom_nodelab(aes(label=label),nudge_x=-.18,nudge_y=.4,size=2)+
  geom_treescale(x=7)
## Warning: Duplicated aesthetics after name standardisation: size
ggpubr::ggarrange(tree_raxml_outgrouproot_visual,tree_compare,nrow=1)

Cladogram, circular tree

tree_compare <- ggtree(tree_raxml_outgrouproot_compare, layout='circular',branch.length = 'none')+
  geom_tiplab(size=2,offset=.01)+
  geom_tippoint()+
  geom_nodelab(aes(label=label),nudge_x=-.18,nudge_y=.4,size=2)+
  geom_treescale(x=20)

ggpubr::ggarrange(tree_raxml_outgrouproot_visual,tree_compare,nrow=1)

Circular tree with branch lengths, all tips aligned to the same branch length

tree_compare <- ggtree(tree_raxml_outgrouproot_compare, layout='circular')+
  geom_tiplab(size=2,align=TRUE,offset=.01)+
  geom_tippoint()+
  geom_nodelab(aes(label=label),nudge_x=-.18,nudge_y=.4,size=2)+
  geom_treescale(x=4)
## Warning: Duplicated aesthetics after name standardisation: size

## Warning: Duplicated aesthetics after name standardisation: size
ggpubr::ggarrange(tree_raxml_outgrouproot_visual,tree_compare,nrow=1)

Exercise 3

  1. Read in two trees: 'raw_enoyl_seqs.nwk' and 'RAxML_bipartitions.raw_enoyl_seqs' with ape::read.tree(). Store in the variables tree1_fasttree and tree2_raxml

  2. For both tree1_fasttree and tree2_raxml, use ape::root() to root the tree to the leaf using outgroup='tr|E6KYQ2|E6KYQ2_9PAST'

  3. Run ggtree() on both both tree1_fasttree and tree2_raxml and store them in variables tree1_fasttree_visual and tree2_raxml_visual

  4. Use ggpubr::ggarrange() to plot both tree1_fasttree and tree2_raxml side by side

  5. Add the annotation layer scale_x_reverse() to flip the orientation of the tree2_raxml_visual tree

  6. Use ggpubr::ggarrange() to plot the flipped tree2_raxml_visual next to tree1_fasttree_visual

One last thing we’ll touch on are zooming in on specific areas of the tree and collapsing branches.

The nodes of each tree can be visualized

tree_compare <- ggtree(tree_raxml_outgrouproot_compare)+
  geom_tiplab(size=2,align=TRUE,offset=.01)+
  geom_nodepoint(color='blue',shape=9)+
  geom_tippoint()+
  geom_nodelab(aes(label=label),nudge_x=-.08,nudge_y=.4,size=2)+
  geom_treescale(x=4)
## Warning: Duplicated aesthetics after name standardisation: size
tree_compare

In each ggtree object there exist a list of node numbers. These numbers are different from the bootstrap numbers assigned to a node.

We can list ALL node numbers (including those that belong to leaves) on the tree using geom_text2()

Importantly, we display the values by using aes(label=node)

ggtree(tree_raxml_outgrouproot_compare)+
  geom_text2(aes(label=node), hjust=-.3,size=3)

This gives us all node numbers. But we really only want those that are at actual nodes at each branching and not those that are on the leaves. To do so, we use geom_nodelab().

Again, we show node values by using aes(label=node)

ggtree(tree_raxml_outgrouproot_compare)+
  geom_nodelab(aes(label=node),hjust=-.3,size=3,color='blue')+
  geom_treescale(x=2)

Knowing these node values is useful because now we can use this information to subset the tree and zoom in on a specific portion of it.

tree_raxml_outgrouproot_visual <- ggtree(tree_raxml_outgrouproot_compare)+
  geom_nodelab(aes(label=node),hjust=-.3,size=3,color='blue')+
  geom_tiplab(size=3)+
  geom_treescale(x=7)

zoomin_visual <- viewClade(tree_raxml_outgrouproot_visual, node=25)

ggpubr::ggarrange(tree_raxml_outgrouproot_visual,zoomin_visual)

If we know the structure of the tree, and its nodes, we can clean up the code so that we show bootstraps instead of internal node numbers.

tree_raxml_outgrouproot_visual <- ggtree(tree_raxml_outgrouproot_compare)+
    geom_nodelab(aes(label=label),nudge_x=-.15,nudge_y=.4,size=2)+
  geom_tiplab(size=3)+
  geom_treescale(x=7)

zoomin_visual <- viewClade(tree_raxml_outgrouproot_visual, node=25)

ggpubr::ggarrange(tree_raxml_outgrouproot_visual,zoomin_visual)

Exercise 4

  1. Read in the 'RAxML_bipartitions.raw_enoyl_seqs' tree using ape::read.tree(). Store in the variable tree1_raxml. Root the tree to 'tr|E6KYQ2|E6KYQ2_9PAST' using ape::root. Make a tree, tree1_raxml_visual with ggtree() showing each internal label=node using geom_nodelab()

  2. Zoom in on node 35 with viewClade() and store the zoomed in visual in the variable zoomin_visual1

Branches can be collapse()’d at specified nodes. The node can be marked with geom_point2() to more clearly indicate the branch was collapsed.

# collapse, no marker
tree_raxml_outgrouproot_visual2 <- collapse(tree_raxml_outgrouproot_visual, node=36)

tree_raxml_outgrouproot_visual2

# collapse, with marker
tree_raxml_outgrouproot_visual2 <- tree_raxml_outgrouproot_visual2+
  geom_point2(aes(subset=(node==36)))

tree_raxml_outgrouproot_visual2

# Multiple nodes collapsed 

tree_raxml_outgrouproot_visual3 <- collapse(tree_raxml_outgrouproot_visual2, node=24)+
  geom_point2(aes(subset=(node==24)))

tree_raxml_outgrouproot_visual3

Nodes can be collapsed but still show the extent of the phylogenetic diversity that they capture

collapse(tree_raxml_outgrouproot_visual, node=36, 'max')

collapse(tree_raxml_outgrouproot_visual, node=36, 'min')

collapse(tree_raxml_outgrouproot_visual, node=36, 'mixed')

Branches can be highlighted using geom_highight()

tree_raxml_outgrouproot_visual4 <- tree_raxml_outgrouproot_visual+
  geom_highlight(35, 'lightblue')

tree_raxml_outgrouproot_visual4

Branches and clades can be annotated using geom_cladelabel()

tree_raxml_outgrouproot_visual+
  geom_cladelabel(node=35, label='fabG', align=TRUE,offset = 2, angle=90, vjust=2, hjust=0.5)
## Warning: Ignoring unknown parameters: vjust

One last thing we might want to do is not display any bootstrap value less than 50. This requires altering the tree itself, not the tree visualization object.

tree_raxml <- ape::read.tree('RAxML_bipartitions.raw_enoyl_seqs')

tree_raxml_outgrouproot <- ape::root(tree_raxml ,outgroup='tr|E6KYQ2|E6KYQ2_9PAST')

tree_raxml_outgrouproot$node.label[as.numeric(tree_raxml_outgrouproot$node.label)<50] <- '' # make bootstrap values less than 50 an empty character
## Warning in
## tree_raxml_outgrouproot$node.label[as.numeric(tree_raxml_outgrouproot$node.label)
## < : NAs introduced by coercion
tree_raxml_outgrouproot$node.label[as.character(tree_raxml_outgrouproot$node.label)=='Root'] <- '' #make the word `Root` an empty character

Now we can combine all of the above steps to make an informative plot.

We know our node numbers from previous visualizations but we want to show bootstrap values. We also want bootstrap values less than 50 to not be displayed since we deem them uninformative and cluttering the visualization.

We also want to highlight that fabG has some variation but that one fabG, which is novel, is on a different clade, closer to FabV. We are not too interested in fabV other than it’s location, so we want to collapse it.

Putting it all together.

tree_raxml_outgrouproot_visual6 <- ggtree(tree_raxml_outgrouproot)+
  geom_tiplab(size=2,offset=0.1)+
  geom_tippoint()+
  geom_nodelab(aes(label=label),nudge_x=-.1,nudge_y=.4,size=2)+
  geom_treescale(x=7, color=NA)

tree_raxml_outgrouproot_visual6 <- collapse(tree_raxml_outgrouproot_visual6, node=24,'max')

tree_raxml_outgrouproot_visual6 <- tree_raxml_outgrouproot_visual6+
  geom_highlight(35, 'lightblue')+
  geom_cladelabel(node=35, label='fabG', align=TRUE,offset = 1.2, angle=90, vjust=2, hjust=0.5)+
  geom_cladelabel(node=24, label='fabV', align=TRUE,offset = -.2, angle=90, vjust=2, hjust=0.5)+
  geom_cladelabel(node=2, label='novel\nfabG', align=TRUE,offset = 1.1, angle=90, vjust=2, hjust=0.5, color='blue')
## Warning: Ignoring unknown parameters: vjust

## Warning: Ignoring unknown parameters: vjust

## Warning: Ignoring unknown parameters: vjust
tree_raxml_outgrouproot_visual6
## Warning: Removed 5 rows containing missing values (geom_point_g_gtree).

There are many more options for manipulating trees available, view them here!

Thanks for coming to Day 1!